Outstanding Questions and To-do's:

1. Usage of x1s2, x1s3->Xs<dir,side>, etc. These are "area-averaged" terms but documentation in Athena++ is not specific what coordinate and over what area the coordinate is averaged. My best guess is that x1s2 is the "1" coordinate averaged over side "2". These terms are used for prolongation of field-centered variables and so might not be needed until a CT MHD method exists in AthenaPK.

2. Naming of h2v, h2f, dh2vd1, dh32vd2, etc. What are these terms? In Athena++, they seem to only be used for viscosity and conduction. Ideally we'd rename these with something more descriptive and replace the different functions with a handful of functions with template parameters like we've done with X<dir> and Xf<dir,face>. However, since their use is limited, we could spin these off into another friend object of the coordinate classes.

3. Do geometric source terms belong within Parthenon? Can these terms be abstracted or are they specific for AthenaPK's implementation?

4. How should coordinate terms be cached to save computation? Both the cylindrical and spherical coordinates in Athena++ make use of cached areas for coordinate systems. Some of these cached arrays may only be used in certain contexts (i.e. just for viscosity). Additionally, the current implementation of cached coordinate data in this PR is likely broken. For now, we could replace the cached data with functions so that for initial testing everything is computed on the fly.

5. Is there a simple test of coordinate systems that we can add to Parthenon? Maybe a disk advection problem?

@forrestglines my apologies---After a longer discussion I think @lroberts36 @jdolence and I have all convinced ourselves that my above comments are not right:

The geometric source terms can be very likely expressed as differences in "areas" or "generalized volumes"

This is still mostly true, I think, but that doesn't mean all geometry output can be dropped. In particular, the face area argument holds for equations similar to the hydro equations, but not for all conservation laws. Equations with a diffusive, e.g., second derivative, term don't work out as nicely.

For prolongation, there are two ways to go. The way we copied from Athena++ uses differences of centroid positions, which requires storing centroids. However, you can also "densitize" the prolongated quantities before and after interpolation. For linear reconstructions, this is equivalent---there aren't that many ways do something physically correct with only 2 degrees of freedom. For higher-order reconstructions, the densitized approach is much cleaner.

I worked out the densitized interpolation vs centroids... I think there's reason to provide centroids---they are not exactly the same as the densitized thing, even in the linear case. My bad on this.

Regarding your questions/TODOs:

Usage of x1s2, x1s3->Xs<dir,side>, etc. These are "area-averaged" terms but documentation in Athena++ is not specific what coordinate and over what area the coordinate is averaged. My best guess is that x1s2 is the "1" coordinate averaged over side "2". These terms are used for prolongation of field-centered variables and so might not be needed until a CT MHD method exists in AthenaPK.

I think in the "topological element" framework that @lroberts36 developed, these become centroid values based on different topological elements, e.g., Face1 or Face2. We probably want to rethink the naming convention here... s for "centroid" and c for cell-center is in my opinion confusing. I would rather switch to X for coordinate location (templated on topological element) and Xc for centroid (templated on topological element).

Naming of h2v, h2f, dh2vd1, dh32vd2, etc. What are these terms? In Athena++, they seem to only be used for viscosity and conduction. Ideally we'd rename these with something more descriptive and replace the different functions with a handful of functions with template parameters like we've done with X

and Xf<dir,face>. However, since their use is limited, we could spin these off into another friend object of the coordinate classes.

I looked a bit more carefully at this... In conduction and viscosity, h*v shows up in calculations of derivatives in combinations with grid spacings that indicate they're used for physical lengths. I.e., the line element in an integral so you can compute the length in a derivative.

I believe these are also expressed more descriptively in other parts of the coordinate machinery. For example, the term

(pco_->h2v(i)*pco_->dx2v(j-1))

that shows up in the x2 flux of conduction in Athena++ can be re-expressed as VolCenter2Length in the Athena++ coordinates class. I wonder if there's a lack of code cleanup here where essentially identical code paths were not unified. It might be worth asking the Athena++ folks.

The dh*vd* versions are the derivatives of these lengths with respect to another direction. And v vs f is for volume (i.e., cell center) vs., cell face. This also shows up in equations with second derivatives, but in particular, it shows up when differentiating vector quantities in those second derivative terms. These are coming from the parallel transport of vectors in curved space. I.e., they represent that a radial vector does not point in the same direction at $\phi = 0$ and $\phi=\pi/2$.

There is a mathematical formalism for describing these relevant quantities. In particular the h* terms are formally described as part of the metric tensor, which describes the physical magnitudes of vectors in curved space and, as a corollary, describes lengths of paths, which are the integrals of the magnitudes of tangent vectors. The dh*vd* type terms are formally described by the Christoffel symbols these are essentially rotation coefficients describing how a vector rotates as you move in a coordinate direction (for example $\theta$).

In full generality the metric is a symmetric 3x3 matrix, so it has six possible nonzero terms. In Cartesian, spherical, and cylindrical coordinates, it is a diagonal matrix, so it's 3 nonzero numbers. In full generality, the Christoffel symbols are three symmetric matrices, so 18 nonzero numbers. However for the Cartesian, spherical, and cylindrical coordinates, there are at most two nonzero terms.

Do geometric source terms belong within Parthenon? Can these terms be abstracted or are they specific for AthenaPK's implementation?

I could go either way on this. I think one option would be to use the above geometric description, using the metric tensor and christoffel symbols. Then to get to the AthenaPK version of those equations, when you get there, it would be straightforward to write a translation layer to represent them as these h and dh type terms. That would be the most general solution from the perspective of a library. Alternatively we could provide them as-is with new names, or not provide them. I lean towards the geometric description. Curious what other people think. (And if folks are interested I could prepare a very brief primer on the fomalism and how it relates specifically to the coordinates stuff in a finite volumes context.)

One thing I wanted to ask about---What is the team's perspective on what parthenon should support broadly? Should parthenon only support Cartesian, Cylindrical, and Spherical, but provide enough hooks for downstream codes to add more generic objects? Or should we design with more generality in mind? Or should we just restrict ourselves to things that are no more complex than Spherical? One reason I ask is that the geometric language I describe above is dramatically simpler for the special cases we are considering here than in full generality. My lenaing is to support only Cartesian, cylindrical, and spherical in parthenon proper, but try to design with letting downstream codes do more general things if they want.

How should coordinate terms be cached to save computation?

This may depend on CPU vs GPU. IMO the reason to cache is that sin and cos are transcendental functions, and thus expensive to compute. But for these simple coordinate systems, we may discover that memory bandwidth vs flops tradeoffs mean that the cache is actually slower.

Additionally, the current implementation of cached coordinate data in this PR is likely broken. For now, we could replace the cached data with functions so that for initial testing everything is computed on the fly.

Agreed. Let's just do on-the-fly and see if we need caching at all.

Is there a simple test of coordinate systems that we can add to Parthenon? Maybe a disk advection problem?

I'd be in favor of that. A more strenuous test might be can we advect a Gaussian blob (with a constant linear velocity) through a spherical/cylindrical coordinate system. I.e., it should still go in a straight line, despite the coordinate system being curved.

For completeness, the metric tensor for cylindrical coordinates is:
$$\gamma = diag(1, r^2, 1)$$
and for spherical coordinates:
$$\gamma = diag(1, r^2, r^2 \sin^2\theta)$$

The h* terms in Athena++ are the so-called scale factors which are the square root of the diagonal components of the metric---they're only defined for diagonal metrics, which is why they're only in the Cartesian, cylindrical, and spherical coordinate objects. I would be supportive of including them over the metric if they are called scale factor rather than h as I think the latter is ambiguous.

The derivatives of the scale factors contain similar information to the christoffel symbols. The Christoffel symbols in cylindrical coordinates have essentially one nonzero term:
$$\Gamma^r_{\theta\theta} = - \Gamma^\theta_{r r} = -\frac{1}{r}$$
In Spherical coordinates, they have the following nonzero terms:
$$\Gamma^r_{\theta\theta} = \Gamma^r_{\phi\phi} = - \Gamma^\theta_{r r} = -\Gamma^\phi_{r\phi} = -\frac{1}{r}$$
and
$$\Gamma^\theta_{\theta\phi} = -\Gamma^\phi_{\theta\theta} = \frac{\cot\theta}{r}$$
I would be open to including the Christoffel symbols themselves or the gradients of the scale factors if they're named descriptively.

Oh one final thought: We might want to add functions that transform coordinate positions to Cartesian positions for the purpose of vis tools... And we may also want Jacobians to transform to Cartesian coordinates for visualizing vectors. Those are just for outputs though... not needed in an actual simulation. And they could be encoded instead in, e.g., yt.

A few comments below. But a bigger picture comment also: would you be willing to chat with myself, Josh, and Luke about what's actually needed here? We were chatting today and we're pretty sure that (at least in a finite volume context) this interface can be simplified.

Thank you so much for looking over this PR. Still early stages and we're open to changes/suggestions.

I'd be open to chatting if we think there's a need. Talking others in T-2 though, our initial guiding principle is "Make it work like Athena++ first". Since we've generally followed Athena++'s AMR strategy and AthenaPK's hydro implementation follows Athena++, following Athena++'s curvilinear strategy hopefully should be the minimal set of changes.

Ideally in the future we'd like swap linear momentum for angular momentum in the conserved variables. I'm not sure how many changes that would require to Parthenon and AthenaPK and might end up in a code separate to AthenaPK.

There is a mathematical formalism for describing these relevant quantities. In particular the h* terms are formally described as part of the metric tensor, which describes the physical magnitudes of vectors in curved space and, as a corollary, describes lengths of paths, which are the integrals of the magnitudes of tangent vectors. The dhvd type terms are formally described by the Christoffel symbols these are essentially rotation coefficients describing how a vector rotates as you move in a coordinate direction (for example ).

The h* terms in Athena++ are the so-called scale factors which are the square root of the diagonal components of the metric

Ah, thank you! This was the insight I was hoping for, seeing as I'm newer to using different coordinate systems.

Since they're only used in specific use cases, it might be better to spin these off into different objects templated on coordinate system i.e. ScaleFactors<class Coords> and Christoffels<class Coords> as we need them. I'm not sure how this would be implemented with cached coordinate data if we decide to do that later (make them member objects of Coordinate classes?) but we can figure that out later.

I think in the "topological element" framework that @lroberts36 developed, these become centroid values based on different topological elements, e.g., Face1 or Face2. We probably want to rethink the naming convention here... s for "centroid" and c for cell-center is in my opinion confusing. I would rather switch to X for coordinate location (templated on topological element) and Xc for centroid (templated on topological element).

Bit of clarification, in this PR Xc becomes "X-centriod" (center of volume) instead of just "X-center" (midpoint of coordinates). For cylindrical the centroid r coordinate differs from the midpoint r and for spherical the centroid r and centroid theta (azimuthal?) differ from their midpoints. As far as I can tell, midpoints don't appear in Athena++, centoid positions are always used. For the topological element framework X function that @lroberts36 made, I think cell-centers should return centoids since we apparently don't use midpoints.

Xs is a "face-averaged" position which is apparently different from the centoid position for only the specific case of the r coordinate in the spherical coordinate system and is only used for prolongation of face-centered fields. I chose Xs because it matches Athena++'s convention of x1s2 and Xfa would confuse it with the "function argument" convention we're already using. I'd be open to suggestions of better names for "face-averaged X"

Bit of clarification, in this PR Xc becomes "X-centriod" (center of volume) instead of just "X-center" (midpoint of coordinates). For cylindrical the centroid r coordinate differs from the midpoint r and for spherical the centroid r and centroid theta (azimuthal?) differ from their midpoints. As far as I can tell, midpoints don't appear in Athena++, centoid positions are always used.

Interesting... I'm not sure, but I think we may want to retain the ability to use coordinate midpoint values in addition to the centroid values. I can see that being valuable in some contexts---for example the volume-weighted interpolation (densitized) if anyone wants to use that. That said maybe that philosophy and the centroid philosophy shouldn't meet.

Curious what other folks think about that.

Xs is a "face-averaged" position which is apparently different than the centoid position for only the specific case of the r coordinate in the spherical coordinate system and is only used for prolongation of face-centered fields.

I'm not sure I understand this---my thinking was that a "face-averaged" position would be just the centroid position based at the face topological element. I was thinking there was sort of a generic generalization of the cell centroids for each topological element which matches that definition. What am I missing here?

I'd be open to chatting if we think there's a need. Talking others in T-2 though, our initial guiding principle is "Make it work like Athena++ first". Since we've generally followed Athena++'s AMR strategy and AthenaPK's hydro implementation follows Athena++, following Athena++'s curvilinear strategy hopefully should be the minimal set of changes.

I think it might be worth touching base at some point. I think your strategy makes sense for AthenaPK. From the parthenon perspective, I think we want the minimal feature set that spans AthenaPK's needs and the needs of other downstream codes that might want to use this machinery.

Ideally in the future we'd like swap linear momentum for angular momentum in the conserved variables. I'm not sure how many changes that would require to Parthenon and AthenaPK and might end up in a code separate to AthenaPK.

Interesting... That probably requires coordinate transformation matrices (i.e., the Jacobian). There may be some need for that anyway in pgens and vis tools though.

Since they're only used in specific use cases, it might be better to spin these off into different objects templated on coordinate system i.e. ScaleFactors and Christoffels as we need them. I'm not sure how this would be implemented with cached coordinate data if we decide to do that later (make them member objects of Coordinate classes?) but we can figure that out later.

That's an interesting idea. Broadly I think this makes sense. We should think about how that fits into the meshblock and pack machinery though to make sure they're accessible where they're needed.

Alternatively, one could imagine caching them on the mesh as variables. Maybe some automated machinery for computing scale factors/christoffels/etc centered at relevant topological elements on the grid? We could just treat them as actual mesh variables? That might cause them to have pretty nonstandard shapes per meshblock though.

I'm currently working on a private branch of AthenaPK using these constructs in Parthenon. Once I have something partially working I can share that private branch to interested LANL employees and make an advection example in Parthenon.

I think we might have enough for cell-centered MHD in this PR already. I'm finding some of the functions in the coordinate class (like the Coord_src terms) might better belong downstream. If we decide to cache these values we'll have to show an example in Parthenon on how to do that.